Moving solitons in the discrete nonlinear Schrodinger equation 
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Using the method of asymptotics beyond all orders, we evaluate the amplitude of radiation from 
a moving small-amplitude soliton in the discrete nonlinear Schrodinger equation. When the non- 
linearity is of the cubic type, this amplitude is shown to be nonzero for all velocities and therefore 
small-amplitude solitons moving without emitting radiation do not exist. In the case of a saturable 
nonlinearity, on the other hand, the radiation is found to be completely suppressed when the soliton 
moves at one of certain isolated 'sliding velocities'. We show that a discrete soliton moving at a 
general speed will experience radiative deceleration until it either stops and remains pinned to the 
lattice, or — in the saturable case — locks, metastably, onto one of the sliding velocities. When the 
soliton's amplitude is small, however, this deceleration is extremely slow; hence, despite losing en- 
ergy to radiation, the discrete soliton may spend an exponentially long time travelling with virtually 
unchanged amplitude and speed. 

PACS numbers: 05.45.Yv, 42.65.Tg, 63.20.Pw, 05.45.Ra 



I. INTRODUCTION 

This paper deals with moving solitons of the discrete 
nonlinear Schrodinger (DNLS) equation. The earliest ap- 
plications of the cubic DNLS equation were to the self- 
trapping of electrons in lattices (the polaron problem) 
and energy transfer in biological chains (Davydov soli- 
tons) - see the reviews [1, 2, 3] for references. Relatedly, 
the equation arises in the description of small-amplitude 
breathers in Frcnkcl-Kontorova chains with weak cou- 
pling [1]. In optics the equation describes light-pulse 
propagation in nonlinear waveguide arrays in the tight- 
binding limit [4, 5]. Most recently the DNLS equation 
has been used to model Bosc-Einstcin condensates in op- 
tically induced lattices [6]. 

The question of the existence of moving solitons in the 
DNLS equation has been the subject of debate for some 
time [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. 
Recently, Gomcz-Gardenes, Fiona, Falo, Peyrard and 
Bishop [20, 21] have demonstrated that the stationary 
motion of pulses in the cubic one-site DNLS (the 'stan- 
dard' DNLS) is only possible over an oscillatory back- 
ground consisting of a superposition of plane waves. This 
result was obtained by numerical continuation of the 
moving Ablowitz-Ladik breather with two commensurate 
time scales. In our present paper, we study the travelling 
discrete solitons analytically and independently of any 
reference models. Consistently with the conclusions of 
[20, 21], we will show that solitons cannot freely move in 
the cubic DNLS equation; they emit radiation, decelerate 
and eventually become pinned by the lattice. We shall 
show, however, that this radiation is exponentially small 
in the soliton's amplitude, so that broad, small-amplitude 
pulses arc highly mobile and are for all practical purposes 



indistinguishable from freely moving solitons. 

In the context of optical waveguide arrays — important 
not only in themselves but also as a first step to un- 
derstanding more complicated optical systems such as 
photonic crystals — interest among experimentalists [22, 
23, 24] has recently shifted away from media with pure- 
Kerr nonlinearity (which gives rise to a cubic term in 
the DNLS equation) and towards photorefractive media, 
which exhibit a saturable nonlinearity [25, 26, 27, 28, 29]. 
In practice such arrays may be optically induced in a pho- 
torefractive material [22, 23] or fabricated permanently - 
see [24], for example. The study of solitons in continuous 
optical systems with saturable noniinearitites has a long 
history; interesting phenomena here include bistability 
[30, 31], fusion [32] and radiation effects [33] which do 
not arise in the cubic equation. As for the discrete case, 
the work of Khare, Rasmussen, Samuclscn and Saxena 
[34] suggests that the saturable one-site DNLS may be 
exceptional in the sense of ref. [35]; that is, despite not 
being a translation invariant system, it supports trans- 
lationally invariant stationary solitons. This property is 
usually seen as a prerequisite for undamped motion in 
discrete equations (see e.g. [36]) and indeed, the numeri- 
cal experiments of Viccncio and Johansson [29] have re- 
vealed that soliton mobility is enhanced in the saturable 
DNLS equation. 

The DNLS equation with a saturable nonlinearity is 
the second object of our analysis here; our conclusions 
will turn out to be in agreement with the numerical ob- 
servations of ref. [29] . We will show that for nonlincar- 
ities which saturate at a low enough intensity, solitons 
can slide — that is, move without radiative deceleration — 
at certain isolated velocities. These 'sliding' solitons are 
examples of embedded solitons. 

The usual saturable DNLS equation is the discrete 
form of the Vinetskii-Kukhtarcv model [37] : 
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In order to encompass both cubic and saturable nonlin- 
earities in a single model, we shall instead consider the 
equation 



i(j) n + 4> n+ \ + <pn-l + 



2\<Pr, 



1 + lA4>n\ 



= 0, 



(2) 



obtained from eq. (1) by making the transformation = 
yflfye^+^^n and letting \i = 2/7. In the form 
above, 1/// represents the saturation threshold of the 
medium [31], which tends to infinity as one approaches 
the pure Kerr (cubic) case of \i = 0. The higher the value 
of /it, the lower is the intensity at which the nonlinearity 
saturates. 

This paper is structured in the following way. In sec- 
tion II, we construct a small-amplitude, broad travelling 
pulse as an asymptotic series in powers of e, its am- 
plitude. The velocity and frequency of this soliton are 
obtained as explicit functions of e and its carrier-wave 
wavenumber. Then in section III, the main section of 
this paper, we derive an expression for the soliton's radi- 
ation tails and measure their amplitude using the method 
of asymptotics beyond all orders. In section IV, we inves- 
tigate the influence of this exponentially weak radiation 
on the soliton's amplitude and speed. Finally, in section 
V, we summarise our work and make comparisons with 
some earlier results. 



where = d n ip/dX n . 

From now on we assume that e is small. Our aim in this 
section is to find an approximate solution to eq. (6) — and 
hence eq. (2) — with ip = 0(e). (That is, we are looking 
for small, broad pulses which modulate a periodic carrier 
wave.) To this end, we expand ip, uj and v as power series 
in e: 



ip = e(ip + eipi H ), 

lu = ojq + e 2 uj 2 + • ■ ■ , 
v = vq + e 2 V2 + • • • . 



(7a) 
(7b) 
(7c) 



(We are not expanding k as we consider it, along with 
e, as one of the two independent parameters character- 
ising our solution.) Substituting these expansions into 
eq. (6) gives a hierarchy of equations to be satisfied at 
each power of e by choosing uj n and v n properly. In non- 
linear oscillations, this perturbation procedure is known 
as Lindstedt's method [38]. 
At the order e 1 we obtain 



ujq = 2 cos k, 



while the order e 2 gives 



2 sinfc. 



(8) 



(9) 



II. ASYMPTOTIC EXPANSION 
A. Leading order 

We begin by seeking solutions of the form 

where 

X = e(n - vt) 



(3) 



(4) 



and e is a parameter. By analogy with the soliton of 
the continuous NLS, we expect the discrete soliton to 
be uniquely characterised by two parameters — e.g. e and 
k — while the other two (to and v) are expected to be 
expressible through e and k. Substituting the Ansatz 
(3), (4) into (2) gives a differential advance-delay equa- 
tion 



ij{X + e)e lk + i>(X - e)e 
— ievip (X) 



— i k 



jiP(X) 



;i>(X) = 0. (5) 



This can be written as an ordinary differential equation 
of an infinite order: 



B=0 



n=0 



ievip' 



LUtp 

2 V> = 0, (6) 



These two relations correspond to the dispersion of linear 
waves. At the power e 3 wc obtain the following nonlinear 
equation for ipQ\ 

cosfc^o ~ ^2^0 + 2\ipo\ Vo = 0. 

This is the stationary form of the NLS equation, which 
has the homoclinic solution 



-00 = aV cos k sech(aX) 



with 



a 2 = ^. 

cos k 

Returning to the original variable ip, we note that the 
amplitude a can always be absorbed into e, the parameter 
in eq. (4) and in front of ipo in (7a). That is, there is no 
loss of generality in setting a = 1 and letting e describe 
the amplitude (and inverse width) of the pulse instead. 
This allows us to set 



u>2 = cos k. 



(10) 



Note that the coefficients in eq. (5) arc periodic func- 
tions of the parameter k with period 2tt; therefore it is 
sufficient to consider k in the interval (— tt, it). Also, (5) 
is invariant with respect to the transformation k — > —k, 
e — > — e, v — > — v; hence it is sufficient to consider positive 
k only. Finally, our perturbative solution does not exist 
if cos k is negative. Thus, from now on we shall assume 
that < k < tt/2. 
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B. Higher orders 



where 



At the order e n+3 , where n > 1, we arrive at the fol- 
lowing equations for the real and imaginary parts of ip n : 



Ci Rei/vi = 
C Im ip n = 



Re/n-ipQ 
cos k 

cos k 



(11a) 
(lib) 



and 



Co 
C l 



-d 2 /dX 2 
-d 2 /dX 2 



2 secli X, 
6 sech 2 X 



[n/2] 
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n — 1 n-1 rt-l-m 

n — rn 

o + E E 2 ^ 

n—m — 



2 sin & 



rn— 1 
n—1 n — l—rn 

+mE e 

m=0 1=0 



m=l 1=0 



E 



2 cos fc 



/ ( 2 i) / 

V n - m -£-2j ~ U2jWn-m-l-2j 



(2j)\ rn - m -^- 2 i 
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2 sin fc 
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The linear nonhomogeneous ordinary differential equa- 
tions (11) must be solved subject to a boundedness con- 
dition. 

The bounded homogeneous solutions of eqs (11a) and 
(lib) (scchXtanhAT and sechX, respectively) corre- 
spond to the translation- and U(l)-invariances of eq. (5). 
Including these zero modes in the full solution of eqs (11) 
would amount just to the translation of ip(X) by a con- 
stant distance in X and its multiplication by a constant 
phase factor. These deformations are trivial, and hence 
we can safely discard the homogeneous solutions at each 
order of e. 

As ^>o's real part is even and its imaginary part is odd 
(zero), i/>i's real and imaginary parts are also even and 
odd, respectively. The same holds, by induction, to all 
orders of the perturbation theory. Indeed, assume that 
■00 7 ipi) ■ ■ ■ j "011-1 bave even real parts and odd imaginary 
parts. Then it is not difficult to verify that the func- 
tion Re/ n _i(X) is even and Im/„_ipf) is odd. Since 
the operators Co and C\ are parity-preserving, and since 
we have excluded the corresponding homogeneous solu- 
tions, this means that ip n has an even real part and an 
odd imaginary part. Finally, the homogeneous solution 
of (lib) being even and that of (11a) being odd, the 
corresponding solvability conditions are satisfied at any 
order. 

Note that since the solvability conditions do not im- 
pose any constraints on v n and uj n , the coefficients v n 



with n > 2 and u) n with n > 4 can be chosen completely 
arbitrarily. 



C. Explicit perturbative solution to order e 

Solving eqs (11) successively, we can obtain the dis- 
crete soliton (7a) to any desired accuracy. Here, we re- 
strict ourselves to corrections up to the cubic power in e. 
The order e 3 is the lowest order at which the saturation 
parameter fj, appears in the solution. On the other hand, 
it is high enough to exemplify and motivate our choice of 
the coefficients in (7b) and (7c). 

Letting n = 1 in eq. (12), we have 

2i 

fo(X) = y smkip'o" - iv 2 ip' . 
The corresponding solution of eq. (11) is 



0i = 



Vcos k 
1 



- sin k sech X tanh X 
2 



V2 sin k I X sech X 

2 V 3 



Here the term proportional to X sech X decays to zero 
as |X| — > oo; however, it becomes greater than ipo(X) 
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for sufficiently large \X\, leading to nonuniformity of the 
expansion (7a). In order to obtain a uniform expansion, 
the term in question should be eliminated. Being free to 
choose the coefficients v n with n > 2, we use this freedom 
to set 

v 2 = - sink. 

This leaves us with 



ipi = — V cos k tan k sech X tanh X. 

After 'distilling' in a similar way the correction ip^, where 
we fix 

U)4 = — cos k 

to eliminate a term proportional to X sechX, we obtain 
ijj = eV cos k | sech X + — e tan k sech X tanh X 
+ — e 2 [4 sech 3 X - 3 sech X 
+ i tan 2 /c (14 sech 3 X - 13sechX) 
+ 4/jcosfc(2sechX-sech 3 X)] + 0(e 3 )}, (13) 



where 



uj = 2 
v = 2 



2! 



1 + 5!'' 



e 1 + 0(e 6 



cos k, 



G(6 4 ) 



sin k. 



(14a) 
(14b) 



D. Velocity and frequency of the discrete soliton 

In the previous subsection we have shown that fixing 
suitably the coefficients u> n and v n can lead to a uniform 
expansion of -0 to order e 3 . Our approach was based on 
solving for ip n explicitly and then setting the coefficients 
in front of the 'secular' terms XsechX to zero. Here, 
we show that the secular terms can be eliminated to all 
orders - and without appealing to explicit solutions. 

Assume that the secular terms have been suppressed 
in all nonhomogencous solutions tp m with mupton-1; 
that is, let 



V> m (A)->C m e x + (e x ) 

as X — > — oo (to = 0, 1, . 



,n-l) 



(15) 



(Since we know that the real part of ip m is even and the 
imaginary part odd, it is sufficient to consider the asymp- 
totic behaviour at one infinity only.) The constant C m 
may happen to be zero for some m in which case the 
decay of ip m will be faster than e x . Our objective is to 
choose uj n and v n in such a way that ij) n (X) will also sat- 
isfy (15). To this end, we consider the function (12). All 



terms which are trilinear in tjjQ, . . . , tpn-i and derivatives 
of these functions decay as e 3X or faster; these terms in 
fn—i cannot give rise to the secular terms proportional 
to XsechX in ip n . On the other hand, the terms making 
up the first line in (12) tend to 



[n/2] 



2 cos k 

(27+2)! -W *" +3 



2j 



[*f!] 
3=1 



2 sin k 



C n -2j + l 



as X — > — oo. These are 'asymptotically resonant' terms 
in the sense that their asymptotics are proportional to 
the asymptotics of the homogeneous solutions of (11a) 
and (lib). It is these terms on the right-hand sides of 
(11a) and (lib) that give rise to the secular terms in the 
solution ip n . The resonant terms will be suppressed if we 
let 



2 cos k 
(27+2)!' 



V 2 j 



2 sin k 
(27+1)!' 



3 > 1. (16) 



After the resonant terms have been eliminated, the 
bounded solution to (11) will have the asymptotic be- 



haviour ip n 



Cr,e x as X 



-oo. By induction, this 



result extends to all n > 0. 

Substituting (16), together with (8)-(10), into eqs (7b) 
and (7c), and summing up the series, we obtain 



2cosfccoshe, 



2 sin k 



sinh i 



(17) 



the frequency and velocity of the discrete soliton param- 
ctcrised in terms of k and e. 

Equations (17) coincide with the expressions [39] of 
the Ablowitz-Ladik soliton's velocity and frequency in 
terms of its amplitude and wavenumber. The difference 
between the two sets of answers is in that eqs (17) per- 
tain to small-amplitude solitons only, whereas Ablowitz 
and Ladik's formulas are valid for arbitrarily large am- 
plitudes. 

Note, also, that the velocity and frequency (17) do not 
depend on the saturation parameter (i. This is in con- 
trast to the stationary (v = 0) soliton of eqs (5) and 
(6) obtained by Khare, Rasmussen, Samuelsen and Sax- 
cna [34]. The soliton of ref. [34] has its amplitude and 
frequency uniquely determined by fi. 



III. TERMS BEYOND ALL ORDERS OF THE 
PERTURBATION THEORY 

A. Dispersion relation for linear waves 

As X — > — oo, the series (7a) reduces to ip(X) = 



y c 



l C n e 



x 



The convergence of the series 



e n+ ijj n (X) for all X would imply, in particular, the 
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convergence of the series e n+1 C n . Therefore, if the se- 
ries (7a) converged, the solution ip(X) would be decaying 
to zero as X — > — oo. However, although we have shown 
that the series e n+1 ip n (X) is asymptotic to all orders, 
it does not have to be convergent. For instance, it is easy 
to see that the series cannot converge for v = [i = and 
any e; since the advance-delay equation (5) is translation 
invariant, this would imply that we have constructed a 
family of stationary solitons with an arbitrary position 
relative to the lattice. This, in turn, would contradict 
the well established fact that the 'standard' cubic DNLS 
solitons can only be centered on a site or midway be- 
tween two adjacent sites [40] (that is, that the 'standard' 
discretisation of the cubic NLS is not exceptional [41]). 
Thus, we expect that the perturbativc solution (7a) sat- 
isfies tp(X) — > as \X\ — > oo only for some special choices 
of u, e and [i. 

Can eqs (5) and (6) have a bounded solution despite 
the divergence of the corresponding series ^ e n+1 C„? To 
gain some insight into this matter, we linearise eq. (5) 
about ip = and find nondecaying solutions of the form 
t[) = e lQX l e where Q is a root of the dispersion relation 



lu = 2cos(fc + Q) + Qv. 



(18) 



It is easy to check that there is at least one such har- 
monic solution [i.e. eq. (18) has at least one root] if v ^ 0. 
These harmonic waves can form a radiation background 
over which the soliton propagates (as suggested by the 
numerics of [20, 21] and the analysis of a similar prob- 
lem for the (j) 4 kinks in [36]). Being nonanalytic in e, 
such backgrounds cannot be captured by any order of 
the perturbation expansion. 

As we will show below, not only the wavenumbers but 
also the amplitudes of the harmonic waves are nonana- 
lytic in e. This phenomenon was first encountered in the 
context of the breather of the continuous </> 4 model, where 
Elconskii, Kulagin, Novozhilova, and Silin [42] suggested 
that the radiation from the breather could be exponen- 
tially weak. Segur and Kruskal [43, 44] then developed 
the method of 'asymptotics beyond all orders' to demon- 
strate that, in the limit e — > 0, such radiation does exist. 
We will use Segur and Kruskal's method to measure the 
magnitude of the radiation background of the travelling 
discrete soliton. 

Qualitatively, the fact that the radiation is not excited 
at any order of the perturbation expansion is explained 
by the fact that the soliton exists on the long length scale 
X, whereas the radiation has the shorter scale X/e. To 
all orders, the two are uncoupled. 

Using the relations (17), we can rewrite (18) as 



cosh e — cos Q 
(sinhe/e)Q — sin Q 



= tan k. 



(19) 



The left hand side is plotted in fig. 1. For e = it has 
minima at multiples of 2ir where the curve is tangent 
to the horizontal axis. For nonzero e, the minima of 
the curve are lifted off the Q axis slightly. The minima 
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FIG. 1: The left hand side of eq. (19) for two values of e. The 
root(s) Q n of the dispersion relation (18) are located where 
this graph is intersected by a horizontal straight line with 
ordinate equal to tan A;. 



with larger values of Q have smaller elevations above the 
horizontal axis; i.e., the minima come closer and closer 
to the Q axis as Q grows. We see from the figure that for 
k > kmlx, where fcmlx ~ 0.22, there is only one radiation 
mode. Note also that the left-hand side of eq. (19) is 
negative for negative Q; since we have assumed that < 
k < 7r/2, this implies that eq. (19) cannot have negative 
roots. 



B. Radiating solitons 

Although we originally constructed the expansion (7a) 
as an asymptotic approximation to a solution which is 
stationary in the frame of reference moving with the ve- 
locity u, it can also represent an approximation to a time- 
dependent solution ip(X,t). Here ip(X,t) is related to 
4> n (t), the discrete variable in eq. (2), by the substitution 
(3): 



<t>n(t) = il>(X,t)e 



ikn-\-iujt 



(20) 



The coefficients %j) n in the asymptotic expansion of 
ip(X,t) will coincide with the coefficients in the expan- 
sion of the stationary solution i/>(X) = ^ e n+1 ipn if the 
time derivatives dt^n he beyond all orders of e and hence 
the time evolution of the free parameters k and e occurs 
on a time scale longer than any power of e _1 . Physically, 
one such solution represents a travelling soliton slowing 
down and attenuating as the Cherenkov radiation left 
in its wake carries momentum and energy away from its 
core. 

Substituting (20) into eq. (2), gives 



iip t + ijj + e lk + ip e lk — u>ip — ievipx + 



2|Vf 



i + mIVI 



> = o, 

(21) 

where ■0 ± = ip(X i e, t). 

We consider two solutions of this equation which 
both have the same asymptotic expansion (7a), denoted 
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ijj s (X,t) and ip u (X,t), such that i/; s (X,t) — > as X — > 
+00 and ip u (X, t) — > as X — ► — 00. Since the difference 
\& = ^ s — -0„ is small (lies beyond all orders of e), and 
since the solution i/> s can be regarded as a perturbation 
of ip u , ^ obeys the linearisation of eq. (21) about ip u to 
a good approximation. That is. 



1 +Ml^«| 2 



- / A' 



cj^ — ive^x 
2 A1 |^| 2 (|V«| 2 ^ + ^**) 



(i + mI^I 2 ) 2 



= 0. 

(22) 



Since tp u = O (e), we can solve eq. (22) to leading order 
in e by ignoring the last two terms in it; the resulting 
solutions are exponentials of the form e l Q x l< L -'^ lt _ \y c 
make a preemptive simplification by setting O = 0. [That 
fl has to be set equal to zero follows from matching these 
exponentials to the far-field asymptotes of the stationary 
'inner' solution; see section HID below. Physically, f2 = 
implies that the travelling pulse will only excite the 
radiation with its own (zero) frequency in the comoving 
frame.] The leading-order solution of eq. (22) is therefore 



5> 



op) 



(23) 



where Q n (n = 1, 2, . . .) are the roots, numbered in order 
from smallest to largest, of the dispersion relation (18). 
(Recall that since we have taken k in the interval [0, tt/2], 

all the roots Q n are positive.) For k > k^l x 0.22, there 
is only one root, Q%. 

Higher-order corrections to the solution (23) can be 
found by substituting the ansatz 

* = + ef[ n \X) + e 2 fi n) {X) + ■ ■ .}e^ x ^ 

n 

+ [M n) W + z 3 9i n) (*) + - ]e-* 2 " X/£ (24) 

n 

into eq. (22), expanding the advance/delay terms 

fi2 (X ^ e ) anc ^ 92 3 (X ^ e ) m Taylor series in e, 
and making use of the asymptotic expansion (7a), (13) 



for ip u . For instance, the first few corrections are found 
to be 



gi n \x) = 



4i cos k 



tanhX, 



2sin(fc + Q n ) — v 

2 cos k 
w + vq n — 2 cos(fc — Q n ) 



(25) 



secli X. 



Since tp u — > as X — > — 00, it follows that ip s — > \t as 
X — > — 00, and hence, once we know the amplitudes A ni 
we know the asymptotic behaviour of tp s . We shall now 
employ the method of asymptotics beyond all orders to 
evaluate these amplitudes. 

C. 'Inner' equations 

Segur and Kruskal's method allows one to measure the 
amplitude of the exponentially small radiation by con- 
tinuing the solution analytically into the complex plane. 
The leading-order term of -0, eVcos/csechX, has singu- 
larities at X = ^ + iim, n = 0, ±1, ±2, .... In the vicin- 
ity of these points, the radiation becomes significant; the 
qualitative explanation for this is that the sech function 
forms a sharp spike with a short length scale near the sin- 
gularity point, and hence there is a strong coupling to the 
radiation modes, unlike on the real axis. The radiation, 
which is exponentially small on the real axis, becomes 
large enough to be measured near the singularities. 

We define a new complex variable y 7 such that ey is 
small in absolute value when X is near the lowest singu- 
larity in the upper half-plane: 



ITT 

ey = X — . 

y 2 



(26) 



The variables y and X are usually referred to as the 
'inner' and 'outer' variables, respectively - the transfor- 
mation to y effectively 'zooms in' on the singularity at 
X = if. We also define u(y) = ip(X) and w{y) = ip*(X). 
Continuing eq. (5) to the complex plane — i.e. substitut- 
ing u(y) for tp(X) and w(y) for tjj*(X) — we obtain, in the 
limit e — > 0, 



e ik u(y + 1) + er ik u{y - 1) - 2cosfcu(y) - 2is\uku'(y) H = 0, (27a) 

1 + /JUTO 

e~ lk w(y + 1) + e lk w(y - 1) - 2coskw(y) + 2ismkw'(y) H = 0. (27b) 

1 + /iwu 

Here we have used the fact that lu — > 2 cos A: and v — > 2sinfc as e — > 0. Equations (27) are our 'inner equations'; they 
arc valid in the 'inner region' —00 < Key < 00, Imy < 0. (The solution cannot be continued up from the real X axis 
past the singularity at y = 0.) 

Solving the system (27) order by order, we can find solutions in the form of a series in powers of y . Alternatively, 
we can make the change of variables (26) in the asymptotic expansion (13) and send e — > 0. This gives, for the first 
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few terms, 



V cos k ^ — - + tan A; — - 

y %y 

= V cos k I — - — tan A; — ^ 

y 2y 



1 7 

-(1 — fi cos k) + — tan 2 k 

1 7 

— (1 — /i cos fc) + — tan 2 fc 



yi 
yd 



(28a) 
(28b) 



We are using hats over u and iu to distinguish the series 
solution (28) from other solutions of eq. (27) that will 
appear in the next section. The asymptotic series (28) 
may or may not converge. We note a symmetry u(—y) = 
—w(y) of the power-series solution. 



D. Exponential expansion 

In order to obtain an expression for the terms which lie 
beyond all orders of y , we substitute (u, w) = (u, w) + 
(6u, Sw) into eqs (27). Since u and w solve the equations 
to all orders in then provided Su and Sw are small, 
they will solve the linearisation of eqs (27) about (u, w) 
for large \y\. 

Formal solutions to the linearised system can be con- 
structed as series in powers of y^ 1 . Because u and w are 
both O (y~ x ), the leading-order expressions for Su and 
Sw as y — ► oo are obtained by substituting zero for u and 
w in the linearised equations. This gives 



Su -> ^2 J n exp(iq„y), 

n 

Sw — ► ^ K n cxp(-iq n y) as y — > oo, 



(29) 



where q n (n = 0, 1, 2, . . .) are the roots of 

cos(fc + q) — cos k + q sin k = 0. (30) 

Note that the roots q n with n > 1 are given by the e — > 
limits of the roots of the dispersion equation (19): q„ = 
lim £ ^ Qn- I n addition, there is a root q a = which does 
not have a Qo counterpart. 

The full solutions (i.e. solutions including corrections 
to all orders in y _1 ) will result if we use the full inverse- 
power series (28) for u and w; these solutions should have 
the form 



oo ,(n) 

5u = Y. Kn Y.-^T expH?„y), 



Sw = K. t 



Cm 



(31a) 

exp(-iq n y). (31b) 



Note that we have excluded the terms proportional to 
e n«,v f rom this ansatz (i.e. set the amplitudes J n to 
zero) as they would become exponentially large on the 
real X axis. [One can readily verify this by making the 
change of variables (26) in cq. (29).] The coefficients , 



(n) 
'2 ' 



(«) 



. and d\ , d. 



, u 2 , • ■ 



. are found recursively when the 
ansatz (31) is substituted into the linearised equations 
and like powers of y^ 1 collected. In particular, the first 
few coefficients are 



» 



2i cos k 



sin(fc + q n ) — sin k ' 
(„) [cos(fc + q n ) — 2 cos k] cos k 
[sin(fc + q n ) — sinfc] 2 



and 



An) 



0. 



An) 



cos k 



cos(fc 



cos k — q n sin k 



(32a) 



(32b) 



where n = 1, 2, 

Having restricted ourselves to considering the lin- 
earised equations for Su and <5u>, we have only taken into 
account the simple harmonics in (29) and (31). Writing 

Su = e 1 Ui + e 2 U 2 H and Sw = e x W x + e 2 W 2 H , 

where e is an auxiliary small parameter (not to be con- 
fused with our 'principal' small parameter e); substitut- 
ing u = u + Su and w = w + Sw in eqs (27), and solv- 
ing order-by-order the resulting hierarchy of nonhomo- 
geneous linear equations, we can recover all nonlinear 
corrections to Su and Sw. The e 2 -corrections will be 
proportional to e~ I ( 9 ™ +9m ) y ; higher-order corrections will 
introduce harmonics with higher combination wavenum- 
bers. Later in this section it will become clear that e is 
actually of the order exp(— nQi/2e) [see eq. (36) below]; 
hence the amplitudes of the combination harmonics will 
be exponentially smaller than that of exp(— iqiy). 

Now we return to the object that is of ultimate interest 
to us in this work - that is, to the function "J of section 
III B representing the radiation of the moving soliton. We 
need to match \P to the corresponding object in the inner 
region. To this end, we recall that $ = ip s — ip Ul where ips 
and ipu are two solutions of the outer equation (5) which 
share the same asymptotic expansion to all orders. In 
the limit e — > 0, the corresponding functions 



u s (y,t) = V> s (X, t) , w s (y,t) = r s (X, t) , 
u u {y,t) = ip u (X,t), w u (y,t) = ipl(X,t) 



(33) 



solve eqs (27) and share the same inverse-power expan- 
sions. We express this fact by writing 



u s (y,t) 
u u {y,t) 



u{y), 



w s (y,t) 
w u (y,t) 



w{y), 

w(y). 
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Therefore, the difference u s — u u (which results from the analytic continuation of the function \I>) can be identified 
with 8u and w s — w u with Sw. Continuing eq. (24) and its complex conjugate gives, as e — > 0, 



and 



i s - u u = V" lim A n (e) exp - 

n x 



We — W. 



u = V lim ^ (e) exp 

*• — J e— *-0 



2e 
2e 



2e 



l + O 



exp{iq n y) 



2 cos k 



1 



uj + vq n - 2 cos(fc - q n ) y 2 



o 4 



exp(-ig„y) (34a) 



l + O 



^limA„(e)exp ( - 



2e 



cxp(-ig„y) 
2 cos k 



lo + vq n - 2 cos(fc - q n ) y 2 



O 4 



exp(ig„y). (34b) 



[Here we have used eq. (25).] Matching (34a) to (31a) 
and (34b) to (31b) yields then K = and 



lim A n (e) exp f - 



= 0, 



lim A* n (e) exp f-^ 1 ) 



(35a) 
(35b) 



for n = 1,2,.... We note that eq. (35a) follows from 
eq. (35b), while the latter equation can be written, sym- 
bolically, as 



A n (e) — ► A'* exp 



2e 



as e 



0. 



(36) 



Our subsequent efforts will focus on the evaluation of the 
constants K n . 

For k greater than k^ix (approximately 0.22), there 
is only one radiation mode and therefore only one pre- 
cxponential factor, K\. For smaller k we note that the 
amplitude of the nth radiation mode, A n , is a factor of 
ex P {^(2n — 2i)} smaller than Ax, the amplitude of 
the first mode. Referring to fig. 1, it is clear that for n > 
3, the difference Q n — Q\ will be no smaller than 7r. As 
for the second mode, it becomes as significant as the first 
one only when k = O (e 2 ) in which case (Q2 — Qi)/ e — 
O(l). But in our asymptotic expansion of section II 
we assumed, implicitly, that k is of order 1 and so the 
case of k = O (e 2 ) is beyond the scope of our current 
analysis. Therefore, for our purposes all the radiation 
modes with n > 2 (when they exist) will have negligible 
amplitudes compared to that of the first mode, provided 
K\ is nonzero and e is small. For this reason, we shall 
only attempt to evaluate K\ in this paper. 



E. Borel summability of the asymptotic series 

Pomeau, Ramani, and Grammaticos [45] have shown 
that the radiation can be measured using the technique 



of Borel summation rather than by solving differential 
equations numerically, as in Segur and Kruskal's original 
approach. The method has been refined by (among oth- 
ers) Grimshaw and Joshi [46, 47] and Tovbis, Tsuchiya, 
Jaffc and Pelinovsky [48, 49, 50, 51], who have applied 
it to difference equations. Most recently it has been ap- 
plied to differential-difference equations in the context of 
moving kinks in 4 models [36]. This is the approach 
that we will be pursuing here. 

Expressing u(y) and w(y) as Laplace transforms 



U(p)e- pv dp 



and 



w(y) = / W(p)e- py dp, 



(37a) 



(37b) 



where 7 is a contour extending from the origin to infinity 
in the complex p plane, the inner equations (27) are cast 
in the form of integral equations 

f(p)U + fJ.[f(p)U] *U *W + U *U *W = 0, (38a) 

f(-p)W + fJ,[f(-p)W] *W *U + W *W *U = 0, 

(38b) 

where 

f(p) = (coshp — 1) cos k — i(sinhp — p) sin k. 
The asterisk * denotes the convolution integral, 



U(p) * W{p) 



Uip-p^Wip^dp! 



where the integration is performed from the origin to the 
point p on the complex plane, along the contour 7. In de- 
riving eqs (38) , we have made use of the convolution the- 
orem for the Laplace transform of the form (37), where 
the integration is over a contour in the complex plane 
rather than a positive real axis. The theorem states that 

u(y)w{y) = ( [U(p) * W(p)]e~ py dp. (39) 
J 7 
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The proof of this theorem is provided in appendix A. 

We choose the contour so that argp — ► ir/2 as \p\ — ► oo 
along 7. In this case we have e~ py — > as \p\ — > 00 for 
all y along any line —00 < Rey < 00 with Imy < 0. 
Therefore, the integrals in (37) converge for all y along 
this line and any bounded U{p) and W(p). 

The function U (p) will have singularities at the points 
where f(p) vanishes while the sum of the double- 
convolution terms in (38a) docs not. Similarly, W(p) 
will have a singularity wherever f(—p) vanishes [while 
the sum of the double-convolution terms in (38b) does 
not]. Therefore, U and W may have singularities at the 
points where 

coshp — 1 = ±i tan k (sinhp — p), (40) 

with the top and bottom signs referring to U{p) and 
W(p), respectively. The imaginary roots of eq. (40) with 
the top sign are at p = —iq n and those of eq. (40) with 
the bottom sign at p = iq n , where q n are the real roots 
of (30) . (We remind the reader that all roots q n are posi- 
tive.) The point p = is not a singularity as both double- 
convolution terms in each line of (38) vanish here. There 
is always at least one pure imaginary root of eq. (40) (and 

only one if k > fcmlx, where kmlx ~ 0.22). 

In addition, there arc infinitely many complex roots. 
The complex singularities of U [complex roots of the top- 
sign equation in (40)] are at the intersections of the curve 
given by 

cosh re L . 9 , / k \2 
q= . , Wl-snTfc — — -cot A; (41a) 
sm k V V smh re J 

with the family of curves described by 

q = k — arcsin f-r-r — sinfc ) + 27m, (41b) 
V sinh K ) 

71=1,2,... , 

and at the intersections of the curve 



COSfiK / . 9 , / K \ 2 , fin \ 

q= : — —\ 1 — sin k I -— — - cot A; (42a) 

sin k V \ sinh re / 

with the family of curves described by 

q = k + arcsin ( -— — sinfc ) + 7r(2n + 1), (42b) 
V sinh re / 

n= -2,-3,-4,... . 

Here re and q are the real and imaginary part of p: 
p = re + iq. The curve (41a) looks like a parabola opened 
upwards, with the vertex at re = q = 0, and the curve 
(42a) like a parabola opened downwards, with the vertex 
at re = 0, q = — 2cotfc. The curves (41b) and (42b), 
on the other hand, look like parabolas for small re but 
then flatten out and approach horizontal straight lines 
as re — > ±00. [In compiling the list of these 'flat' curves 
in (41b) and (42b), we have taken into account that the 
curve (41b) with n = does not have any intersections 



with the parabola (41a), and the curve (42b) with n = — 1 
does not have any intersections with the parabola (42a) .] 
As k is reduced, the vertex of the parabola (42a) moves 
down along the q axis; the intersections of this parabola 
with the 'flat' curves (42b) approach, pairwise, the q axis. 
After colliding on the q axis, pairs of complex roots move 
away from each other along it. [The parabola (41a) does 
not move as k is reduced, but only steepens, which results 
in the singularities in the upper half-plane approaching 
the imaginary axis but not reaching it until k = 0.] In 
a similar way, the complex singularities of W(p) move 
onto the imaginary axis as k is decreased. In section 
HIE below we will use the fact that the distance from 
any complex singularity to the origin is larger than 27r; 
this follows from the observation that the closest points 
of the curves (41b) and (42b) to the origin are their in- 
tersections with the q axis. These are further away than 
2tt from the origin. 

In addition to singularities at p = —iq n , the func- 
tion U{p) will have singularities at points p = iq n , 
n = 1,2,.... These are induced by the cubic terms in 
eq. (38a); for instance, the singularity at p = iq\ arises 
from the convolution of the term proportional to p in 
U * U with the function W which has a singularity at 
p = iqi. [That U(p) has singularities at p — iq n can 
also be seen directly from eq. (31a).] Similarly, the func- 
tion W(p) will have singularities at points p = —iq n , 
n = 1,2,.... By virtue of the nonlinear terms there 
will also be singularities at the 'combination points' 

i(±q n ± 9m), i(±q n ± 9m ± qj), etc. 

The formal inverse-power series (28), which we can 
write as 

u(y) = ^2pe-ffT' W = X>*7ZTr> (43) 
£=0 y e.=o y 

result from the Laplace transformation of power series 
for U(p) and W(p): 

00 00 

u( P ) = J2 prf> w (p) = E v 'p e - ( 44 ) 

£=0 e=o 

The series (44) converge in the disk of radius qi , centred 
at the origin, and hence can be integrated term by term 
only over the portion of the contour 7 which lies within 
that disk. However, by Watson's lemma, the remaining 
part of the contour makes an exponentially small contri- 
bution to the integral and the resulting series (43) are 
asymptotic as y — > 00. The functions u(y) and w(y) de- 
fined by eqs (37) give the Borel sums of the series u(y) 
and w(y). 

Consider now some horizontal line in the inner re- 
gion; that is, let Imy < be fixed and Rey vary from 
—00 to 00. If the integration contour 7 is chosen to lie 
in the first quadrant of the complex p plane, the func- 
tions u{y) and w(y) generated by eqs (37) will tend to 
zero as Re y — > +00 along this line. Similarly, if it is 
chosen to lie in the second quadrant, they will tend to 
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FIG. 2: The integration contours 7 S and 7„ used to generate 
the solutions (u s ,w a ) and (u u ,w u ) respectively via eqs (37). 
The dots are singularities of W(p). Shown is the situation 
where the linear dispersion relation (30) has only one real 
root, qi. 



zero as Key — > — oo. Suppose there were no singulari- 
ties between two such contours: then the one could be 
continously deformed to the other without any singular- 
ity crossings; i.e. they would generate the same solution 
which, therefore, would decay to zero at both infinities. 
[That is, the oscillatory tails in (31) would have zero am- 
plitudes, K n = 0.] In general, however, U{p) and W{p) 
have singularities both on and away from the imaginary 
axis. In order to minimise the number of singularities 
to be crossed in the deformation of one contour to the 
other, we choose the contours to lie above all singular- 
ities with nonzero real part. (That this is possible, is 
shown in appendix B.) Note also that the imaginary 
part of the singularity grows faster than its real part and 
hence argp should tend to ir/2 as \p\ — > oo along 7; this 
was precisely our choice for the direction of the contours 
7 in the beginning of section III E. 

Let 7,5 and 7 U be two contours chosen in this way, 
with 7 S lying in the first quadrant and 7„ in the second 
quadrant. The solutions u(y), w(y) generated by eqs (37) 
with the contour 7,,. will tend to zero as Re y — > 00 (with 
Imy < fixed). Hence they can be identified with so- 
lutions u s ,w s obtained by the continuation of the outer 
solution tp s which has the same asymptotic behaviour. 
Similarly, the solutions generated by eqs (37) with the 
contour 7„ coincide with solutions u u , w u - like u u , 
w u , the solutions generated by eqs (37) tend to zero as 
Key — > —00 (with fixed Imy < 0). 

Consider, first, solutions w s and w u . Since the con- 
tours 7 S and 7„ are separated by singularities of W(p) 
on the positive imaginary axis, they cannot be continu- 



ously deformed to each other without singularity cross- 
ings and so the solution w s does not coincide with uu u , 
unless the residue at the singularity happens to be zero. 
If we deform 7 S to 7 S and 7„ to "f' u as shown in fig. 2, 
without crossing any singularities, then the only differ- 
ence between the two contours is that 7^ encircles the 
singularities, whereas j' u does not. Therefore, the dif- 
ference w s — w u can be deduced exclusively from the 
leading-order behaviour of W(p) near its singularities. 
There can be two contributions to this difference: the 
first arises from integrating around the poles and is a 
sum of residues, while the second arises if the singularity 
is a branch point. 

To find the singularity structure of the function W(p), 
we equate 

w s -w u = f W(p)e- py dp - I W{p)e- py dp (45) 

to the expansion (31b). The first term in (31b), e~ lqnV , 
arises from the integration of a term (27ri) _1 (p— i(j„) _1 in 
W(p). For such a term, the difference of the two integrals 
in (45) reduces to an integral around a circle centered on 
the point p = iq n : 



1 



27ri / p — iq 



e vv dp = res ■' 



P - iq n 



The term y m e lqVn in (31b), with m = 1, 2, . . ., arises 
from the integration of a term 



1 (p-iq n ) m - 
2iri (m-lY. 



ln(p - iq n ) 



in W(p). This time, p = iq n is a branch point. After 
going around this point along the circular part of "f' s , the 
logarithm increases by 27rz and the difference between the 
integrals in (45) is given by 



{p-iq n ) m - X eT m dp 



[m-iy.Jc 



(m- 1)1 Jo 



z m - 1 e-* v dz. 



where C is the part of 7^ extending from p = iq n to 
infinity. This equals exactly y- m e~ lqnV . 

Thus, in order to generate the full scries (31b) we must 
have 



W(p) 



1 



E 



» 



P - iq n ^-1 {m - 1) 

m— 1 v / 



x (p - iq n ) m 1 Ln(p - iq„) 



W res (p), (46) 



where W reg denotes the part of W which is regular at 
1,2,.. .. By the same process, matching 



P = tq r 



reg 
/) 
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u u to Su in (31a) yields 



oc 

^) = ^E*«E 



27Ti ^ (m- 1)! 

n m=l 

x (p - ife)" 1 " 1 In(p - ig„) + U ieg (p). (47) 

The solution to eqs (38) is nonuniquc; for in- 
stance, if {U(p),W(p}} is a solution, then so is 
{eP*°+<°U(p), e pyo ~ c ° W(p)} with any complex y and Co- 
Also, if {U(p),W{p)} is a solution, {W(-p),U(-p)} is 
another one. We will impose the constraint 



U(p) = W(-p); 



(48) 



this constraint is obviously compatible with eqs (38). 
It is not difficult to see that the reduction (48) singles 
out a unique solution of eqs (38). The motivation for 
imposing the constraint (48) comes from the symmetry 
u{—y) = —w(y) of the power-series solution of eq. (27). 
Using this symmetry in eq. (43), we get pi = (— 1) vt and 
then eq. (44) implies (48) . 

In view of (48), the singularities of U(p) in the up- 
per half-plane are singularities of W(p) in the lower half- 
plane, which fall within W re g(p)t and vice versa. Thus 
we have, finally, 



2ni p — iq n 2iri 



' (to - 1)! 



- (-l) ro dfc>(p + iqn)™- 1 \n{p + iq n )} + W leg (p), (49) 

where W leg (p) is regular at p = ±iq n . We also mention 
an equivalent representation for (49) which turns out to 
be computationally advantageous: 



^ oo 



m=0 



» 



(-l) m d 



(») 



p - iq n p + iq„ 

+ W Tes (p). (50) 



Here D 1 is an integral map: 

D-'fip)^ [ P f( Pl )d Pl ; 

the notation D~ m f(p) should be understood as 

rP rPi rP2 fPm-i 
D~ m f(p) = I dpx / dp 2 / dp a --- dp m f{ Pm ). 



We have also introduced c. 



00 



1 and do 



for econ- 



omy of notation. The only difference between eqs (49) 
and (50) is that the double-sum term on the right-hand 
side of (50) includes some terms which are regular at 
p = ±iq n , whereas in eq. (49), all regular terms are con- 
tained in W rcg (p). 



The residues K n at the poles of W{p) are known as the 
Stokes constants. The leading-order Stokes constant K\ 
can be related to the behaviour of the coefficients in the 
power-series expansion of W(p). Indeed, the coefficients 
in the power series (44) satisfy 



lYdtt (£■ 



to)! 



2nq 1 (iq i y- 



fl 



as £ 



oo. 



(51) 

This is obtained by expanding the singular part of the 
expression (50) in powers of p. (Coefficients of the regular 
part become negligible in the limit £ — > oo compared to 
those of the singular part.) Note that we have ignored 
singularities with nonzero real part and singularities on 
the imaginary axis other than at p = ±iqi- The reason 
is that all these singularities are further away from the 
origin than the points Hqi (in particular all complex 
singularities are separated from the origin by a distance 
greater than 2n), and their contribution to vg becomes 
vanishingly small as £ — ► oo. We have also neglected 
singularities at the 'combination points' because of their 
exponentially small residues. The coefficients vi can be 
calculated numerically; once they are known, it follows 
from (51) that 

R\ = 2ir qi lim [iqx) l v t (52) 



(where we have recalled that c^ 1 ' = 1 and d^ = 0). 

We now turn to the numerical calculation of the coef- 
ficients 



F. Recurrence relation 

To make our forthcoming numerical procedure more 
robust, we normalise the coefficients in the power series 
(44) by writing 



vi = -i 



(53) 



Substituting the expansions (44) with eq. (53) as well 
as the constraint (48) into either of equations (38) and 
equating coefficients of like powers of p, yields the follow- 
ing recurrence relation for the numbers S n (n > 0): 



E 

m=0 



q?6 n . 



(m + 2) 



■n ri 



E 



$n. — r 



J=2 



(n + 2)(n + 1) 



3i< 



n—m—j 



j 




u m—j Vj 



j\(m — j)\ \ m\(n — to)! 



TO! 



(54) 



Here 



(— l)™cosfc, for to even, 
(— l) -1 ^— sin fc, for m odd. 
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Solving eq. (54) with n = gives 5o = \/l — c 2 . There- 
after it can be solved for each member of the sequence 
{Sn} in terms of the preceding ones, and thus each S n 
can be calculated in turn. [Since all the coefficients in 
the recurrence relation (54) are real, the sequence {5 n } 
turns out to be a sequence of real numbers.] 

Once the sequence {S n } has been generated, expression 
(52) can be used to calculate the Stokes constant K\ : 



K x = -2iriqiSe, 



(55) 



for sufficiently large £. Unfortunately, the convergence of 
the sequence {S n } is slow and thus the above procedure 
is computationally expensive. The convergence can be 
accelerated by expanding eq. (51) in powers of small 1/t: 



27T<Jl 



whence 



c 



& 

C- 



-250 



-260 



-270 



-280 




1000 



FIG. 3: Convergence of the sequence on the right-hand side 
of (55) (dashed line) and the 'accelerated' sequence defined 
by the right-hand side of (57) (solid line). Shown are the £th 
approximations to the Stokes constant K\ [the fth members 
of the sequences (55) and (57)] divided by i to get a real value. 
In this plot, n — and k = 0.5. 



c (1) 

-Iviq^i 1 - 

ii 



r- 



o 



(56) 



According to eq. (56), eq. (55) gives K\ with a relatively 
large error of order 1 /£. On the other hand, a two-term 
approximation 



K\ = —2iriqi5e 



(57) 



is correct to 0(l/£ 2 ). More precisely, the relative error 
associated with the answer (57) is given by 



2 L 2 



(4 1) ) 2 + (-i) f 4 1) 



(58) 



In our calculations, we set £ / K\ = 10 5 . Since Ci \ < 

7(1) 



(1) 
2 

and di^' are known constants [given by (32a) and (32b)], 
eq. (58) tells us what I we should take - i.e. how many 
members of the sequence {S n } we should calculate in or- 
der to achieve the set accuracy. Figure 3 illustrates the 
convergence of the approximate values of K\ calculated 
using (55) and (57) as £ is increased. Note the drastic 
acceleration of convergence in the latter case. 

Figure 4(a) shows the calculated Stokes constant as a 
function of k for various values of the saturation param- 
eter /i. First of all, Ki(k) does not have any zeros in 
the case of the cubic nonlincarity (/i = 0). This means 
that solitons of the cubic one-site discrete NLS equation 
[eq. (2) with \i = 0] cannot propagate without losing en- 
ergy to radiation. For /i = 3 the Stokes constant does 
have a zero, but at a value of k smaller than fcmax (where 
fcmax ~ 0.22 = 0.077r). Since higher radiation modes do 
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FIG. 4: The Stokes constant Ki for various values of fi. Note 
the logarithmic scale on the vertical axis. The downward 
spikes extend all the way to — oo; hence each spike corresponds 
to a zero crossing. Each panel shows only a portion of the full 
range < k < n/2; there are no additional zero crossings in 
the part which is not shown. 
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exist in this range of k, there will still be radiation from 
the soliton - unless the 'higher' Stokes constants K 2 (k), 
Ks{k),. . . , happen to be zero at the same value of k. Fi- 
nally, for fi = 4 the zero is seen to have moved just above 
fcmax and for /it = 6 it has an even higher value. There 
are no Q2, Q3, ■ ■ ■ radiations for these k; hence the zeros 
of Ki(k) define the carrier wavenumbers at which the 
soliton 'slides' - i.e. travels without emitting any radia- 
tion. Equation (17) then gives the corresponding sliding 
velocities, for each e. 

Figure 4(b) shows the Stokes constant Ki(k) for higher 
values of the parameter \i. For \i = 12 a second zero of 
the Stokes constant has appeared while for \i = 25, the 
function Ki(k) has three zeros. As /1 is increased, the 
existing zeros move to larger values of k while new ones 
emerge at the origin of the k axis. 

G. Radiation waves 



the first radiation mode, extending to —00, carries energy 
away from the soliton. 

The even-numbered radiation modes (where present) 
in our asymptotic solution (59) have positive group ve- 
locities and hence describe the flux of energy fed into 
the system at the left infinity. A more interesting situa- 
tion is obviously the one with no incoming radiation; the 
corresponding solution is obtained by subtracting off the 
required multiple of the solution of the linearised equa- 
tion, e.g. e lQ2X / e . One would then have a pulse leaving 
the odd modes in its wake and sending even modes ahead 
of it. 

If the first Stokes constant Ki(k) has a zero at some 
k = ki while Q\ is the only radiation mode available (as 
happens in our saturable model with [i greater than ap- 
proximately 4), then according to eq. (59), the radiation 
from the soliton with the carrier-wave wavenumber k\ is 
suppressed completely. 



For not very large \X\, the solution ip s is close to the 
localised pulse found by means of the perturbation ex- 
pansion in section II. As X — ► +00, it tends to zero, 
by definition, while the X — > —00 asymptotic behaviour 
is found from ip s = ip u + iff . Here the solution ip u de- 
cays to zero as X — > —00 and hence tp s approaches the 
oscillatory waveform iff given by eqs (24), (25), and (36): 



mx)^J2 k: > 



-7rQ„/2e 



1 



Ai cos k tan X 
2sin(A; + Q n ) — v 



Oit? 



,iQ n X/t 



(59) 



as X — > —00. Equation (59) describes a radiation back- 
ground over which the soliton is superimposed. As we 
have explained, we can ignore all but the first term in 
the sum. 

To determine whether the radiation is emitted by the 
soliton or being fed into it from outside sources, we con- 
sider a harmonic solution ip = e }Q x h- l ^ lt of the linearised 
eq. (21); the corresponding dispersion relation is 

fi(G) = -2 cos(fc + Q) + lu - Qv. 

The radiation background iff consists of harmonics with 
51 = and Q = Q n where Q n are roots of eq. (18). The 
group velocities of these harmonic waves are given by 



fi'(fi„) = 2sin(fc + Q„) - v. 



(60) 



The <2„'s are zeros of the function f2(Q) and the group 
velocities are the slopes of this function at its zeros; thus 
the group velocities f2 / (Qi), Qf{Q-z)y ■ ■ , have alternating 
signs. The first one, which is the only one that concerns 
us in this work, must be negative. Indeed, the value 

O(0) = — 2cosfc + uj = 2cosfc(coshe — 1) 

is positive, and hence the slope of the function fl(Q) as 
it crosses the Q axis at Qi > is negative. Therefore, 



IV. TIME EVOLUTION OF A RADIATING 
SOLITON 

A. Amplitude-wavenumber dynamical system 

To find the radiation-induced evolution of the travel- 
ling soliton, we use conserved quantities of the advance- 
delay equation associated with eq. (2). In the reference 
frame moving at the soliton velocity v this equation reads 



iip t + f{x+l, t)+(p(x — f, t) — iv(p x + 



2MV 



1 



0. (61) 



The discrete variable 4> n (t) in eq. (2) is related to the 
value of the continuous variable (p(x,t) at the point 
x = n — vt: 4> n (t) = (p(n — vt,t). For future use, we 
also mention the relation between tp(x, t) and the corre- 
sponding solution of eq. (21): 

<p{x, t) = ip(X, t)e ik ( x+v V+ iu)t . (62) 

We first consider the number of particles integral: 



N = 



\ip\ 2 dx. 



Multiplying eq. (61) by tp*, subtracting the complex con- 
jugate and integrating yields the rate of change of the 
integral N: 



dN 
'~dt 



(<p + <p* — c.c.)dx 



a-l 



6+1 



(ip (p* — c.c.)dx + iv\tp\' 



(63) 



In eq. (63), ip = (p(x ± 1, t) and c.c. stands for the com- 
plex conjugate of the immediately preceding term. The 
integration limits a < and b > arc assumed to be 
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large (\a\,b ^> e 1 ) but finite; for example one can take 
a,b = 0{e- 2 ). 

Consider the soliton moving with a positive velocity 
and leaving radiation in its wake. This configuration is 
described by the solution ip s of eq. (21); the correspond- 
ing solution of eq. (61) has the asymptotic behaviour 
ip s — > as x — > +oo. Substituting the leading-order ex- 
pression (23) for the soliton's radiation tail into eq. (62) 
yields the asymptotic behaviour at the other infinity: 



Taking time derivatives of N and P above and discard- 
ing e 1 corrections to e° terms, we deduce that 



iV + tanfc(P-fciV) 

2 cos k 
P-kN 



2e cos k 



tp s (x,t) 



A n [l + 0(e) 



(64) 



as x — > — oo. Since, as we have explained above, the 
Qi radiation is dominant, it is sufficient to keep only the 
n = 1 term in (64). 

Substituting (64) into (63) and evaluating the integral 
over the region (a — 1, a) in the soliton's wake, we obtain 



^L = \K 1 \ 2 e~^^l2sm(k+Q 1 )-v], 



(65) 



where we have used eq. (36). Note that [2 sin(A;+ Qi) — v] 
is the group velocity of the first radiation mode, fi'(Qi), 
which, as we have established, is negative; hence N < 0. 
We now turn to the momentum integral, 



P 



((p*tp - tp x (p*)dx. 



Multiplying (61) by ip*, adding its complex conjugate 
and integrating, gives the rate equation 



dP 
dt 



— = / (tp + <p* x + c.c.)dx 



a-l 



6+1 



(ip tp* x + c.c.)dx 



a-l 



6+1 



-M 2 -^ln(l 



(66) 



Evaluating the right-hand side of (66) similarly to the 
way we obtained (65) and substituting eq. (18) for cj, pro- 
duces 

rlP 

^- = [K^e-^/^k + Qi)[2sin(fe + Q x ) - v}. (67) 

Using the leading-order term of the perturbative solu- 
tion (13) and eq. (62), we can express N and P via e and 
k: 



N = 2e cos k - 
P = 2ke cos k 



In calculating N and P we had to integrate from X = ea 
to X = eb. Since the integrands decay exponentially, and 
since a < and b > were assumed to be much larger 
than e _1 in absolute value, it was legitimate to replace 
these limits with — oo and oo, respectively. The error 
introduced in this way is exponentially small in e. 



Finally, substituting for N and P from cqs (65) and (67), 
we arrive at the dynamical system 



e = \Kr{k)fe-^l< n'(&) I+^I^ (68a) 



k=\K 1 (k)\ 2 e-^l* n'(Q!) 



2 cos k 

Si 



2e cos k ' 



(68b) 



where the group velocity ^l'(Qi) = 2 sin(fc + Q{) — v with 
v = 2(sinhe/e) sinfc, and Qi = Qi(e,k) is the smallest 
root of eq. (19). 
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FIG. 5: The phase portrait of the system (68) in the case 

(a) where the Stokes constant K\ (k) does not have zeros, and 

(b) where Ki(k) has one zero. In (b), the dashed line is the 
line of nonisolated fixed points k = ki. Note that in (b), the 
phase portrait has been replotted on the (e, v) plane; hence 
the dashed line gives the value of the sliding velocity for each 
value of e. In (a), (j, = 0; in (b), fi = 6. 
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B. Soliton's deceleration and sliding velocities 

The vector field (68) is defined for k > , where fc^v 
is the value of k for which the roots Q\ and Qi merge in 
fig. 1 (i.e. the smallest value of k for which the roots Q\ 
and 0,2 still exist). When k = k^ n , the group velocity 
^'(Qi) becomes equal to zero; therefore, the equation 

k = ^^(e) defines a line of (nonisolated) fixed points 
of the system (68). Assume, first, that the saturation 
parameter /i is such that Ki(k) does not vanish for any 

k. For k > fc^in, the factor f2'(Qi) in (68b) is negative; 
since Qi and cosfc are positive for all < k < n/2, 
the derivative k satisfies k < for all times. Hence, all 

trajectories are flowing towards the line k = fc£i n (e) from 
above. It is worth noting that the derivative e is also 
negative, and that the k axis is also a line of nonisolated 
fixed points. However, no trajectories will end there as 
follows from the equation 

de T 1 

Tk= e r U fc+ &(M. ' 

which is a consequence of (68a) and (68b). 

Representative trajectories are plotted in fig. 5(a). 
Since k^ n (e) f=a e 2 /47r is very small for small e, the line 

k = k^J n (e) is practically indistinguishable from the hor- 
izontal axis. Therefore, we can assert that, to a good 
accuracy, k — > as t — > oo. This means that the soliton 
stops moving - and stops decaying at the same time. 

For k smaller than k^- m , the vector field (68) is un- 
defined and we cannot use it to find out what happens 
to the soliton after k has reached fc^ in . The reason for 
this is that the soliton stops radiating at the wavenumbcr 
Qi as k drops below k^ in . In fact our analysis becomes 
invalid as soon as k becomes C(e 2 ) — i.e. even before k 
reaches — because we can no longer disregard the 
n = 2 radiation here. At the qualitative level it is ob- 
vious, however, that the parameter k should continue 
to decay all the way to zero, in a cascade way. First, 
the n = 2 radiation will become as intense as the n = 1 
mode when fc approaches e 2 /47r. Subsequently i.e. for k 
smaller than e 2 /4-7r — the n = 3 harmonic will replace the 
radiation with the wavenumbers Q± and Q2 as a domi- 
nant mode. The n = 4 mode will become equally intense 
near k = fc^ n e 2 /87r; as k drops below e 2 /8ir, both Q3 
and Q4 will cede to Q5, and so on. 

If /.t is such that the Stokes constant K\{k) vanishes 
at one or more values of k, the system has one or more 
lines of nonisolated fixed points, k ~k{. The correspond- 
ing values of v, v = i>i(e) = 2(sinhe/e) sinfe^, define the 
sliding velocities of the soliton - i.e. velocities at which 
the soliton moves without radiative friction. One such 
velocity is shown by the dashed line in fig. 5(b). The 
fixed points (e, k{) are semistable; for k above ki 1 the flow 
is towards the line k = k t but when k is below fcj, the 
flow is directed away from this straight line [sec fig. 5(b)]. 
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FIG. 6: The soliton's decay rate as a function of its velocity, 
for fixed e = 0.1. (Here fj. = 0.) In the inset, the same curve 
is replotted using a logarithmic scale on the vertical axis. 



Since both e and k are negative, the soliton's velocity 
v = 2(sinh e/e) sin k will generally be decreasing - until 
it hits the nearest underlying sliding velocity Vi and locks 
on to it. The ensuing sliding motion will be unstable; a 
small perturbation will be sufficient to take the soliton 
out of the sliding regime after which it will resume its 
radiative deceleration. However, since k is proportional 
to the square of the Stokes constant and not to K\ (k) it- 
self, small perturbations Sk will be growing linearly, not 
exponentially in t. As a result, the soliton may spend 
a fairly long time sliding at the velocity Vi. It is there- 
fore not unreasonable to classify this sliding motion as 
metastable. 

We conclude that the soliton becomes pinned (i.e. 
k — > and so v — > 0) before it has decayed fully (i.e. 
before the amplitude e has decreased to zero). The ex- 
ponential dependence of e and k 011 1/e implies that tall, 
narrow pulses will be pinned very quickly while short, 
broad ones will travel for a very long time before they 
have slowed appreciably. Next, if \x is such that there is 
one or more sliding velocity available in the system, and 
if the soliton is initially moving faster than some of these, 
its deceleration will be interrupted by long periods of un- 
damped motion at the corresponding sliding velocities. 

It is worth reemphasising here that if the amplitude e 
is small, then even if the soliton is not sliding, its decel- 
eration will be so slow that it will spend an exponentially 
long time travelling with virtually unchanged amplitude 
and speed. The deceleration rate —v/v is shown in fig. 6, 
as a function of the soliton's velocity v, for fixed ampli- 
tude e. Note that the decay rate drops, exponentially, as 
the velocity is decreased; this drop is due to the expo- 
nential factor e _7rSl / c in eq. (68). As the velocity v (and 
hence the wavenumber k) decreases, the root Qi(e, k) 
grows towards the limit value of approximately 2-7T (see 
fig. 1). This variation in Qi is amplified by the division 
by small e and exponentiation in e _lr ^ 1 / c . 
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V. CONCLUDING REMARKS 
A. Summary 

In this paper, we have constructed the moving discrete 
soliton of the saturable NLS equation (2) [and hence (1)] 
as an asymptotic expansion in powers of its amplitude. 
The saturable nonlincarity includes the cubic NLS as a 
particular case (for which \x = 0). Our perturbation pro- 
cedure is a variant of the Lindstedt-Poincarc technique 
where corrections to the parameters of the solution are 
calculated along with the calculation of the solution it- 
self. Although the resulting asymptotic series (13) for 
the soliton is generally not convergent, the associated ex- 
pansions for its frequency and velocity sum up to exact 
explicit expressions (17). 

From the divergence of the asymptotic series (13) it 
follows, in particular, that the travelling discrete soli- 
ton does not decay to zero at least at one of the two 
infinities. Instead, the soliton approaches, as X — > oo 
or X — > — oo, an oscillatory resonant background (23) 
where the amplitudes A n of its constituent harmonic 
waves lie beyond all orders of e. For the soliton moving 
with a positive velocity and approaching zero as X — > oo, 
the oscillatory background at the left infinity represents 
the Cherenkov radiation left in the soliton's wake. To 
evaluate the amplitudes of the harmonic waves arising 
as X —>■ — oo, we have continued the radiating soliton 
into the complex plane, where it exhibits singularities. 
We then matched the asymptotic expansion (34) of the 
background near the lowest singularity on the imaginary 
X axis to the far-field asymptotic expansion (31) of the 
background solution of the 'inner' equation - i.e. of the 
advance-delay equation 'zoomed in' on this singularity. 
The asymptotic expansions here are in inverse powers 
of the zoomed variable, y. The amplitudes of the ra- 
diation waves were found to be exponentially small in 
e, with the pre-exponential factors (the so-called Stokes 
constants) being dependent only on the soliton's carrier- 
wave wavenumber. Representing solutions to the inner 
equation as Borel sums of their asymptotic expansions, 
the Stokes constants can be related to the expansion coef- 
ficients; we have calculated these coefficients numerically, 
using algebraic recurrence relations. 

The upshot of the calculation of the leading Stokes 
constant K\ is that in the case of the cubic nonlinearity 
(i.e. for fj, = 0), K\(k) does not vanish for any k. This 
means that the cubic discrete soliton cannot 'slide' - i.e. 
cannot move without radiative friction. The saturable 
solitons, on the other hand, can slide provided the satu- 
ration parameter \i is large enough. This is because, for 
a sufficiently large /j,, the Stokes constant K\{k) is found 
to have one or more zeros k\,hi,.... Since the soliton 
with wavenumber k > 0.22 can have no higher-order res- 
onances, those zeros which satisfy fcj > 0.22 do define the 
wavenumbers at which sliding motion occurs. For each 
value of the soliton's amplitude e, the formula (17) then 
gives the sliding velocities v(fc,,e). 



The calculation of asymptotics beyond all orders is use- 
ful not only for determining the sliding velocities. Know- 
ing the radiation amplitudes has allowed us to derive a 
two-dimensional dynamical system (68) for the soliton's 
parameters. Trajectories of this dynamical system de- 
scribe the evolution of the soliton travelling at a generic 
speed. The evolution turns out to be simple: If /j, is such 
that the Stokes constant K\(k) docs not vanish anywhere 
in the region k > 0.22, the soliton decelerates, although 
very slowly. Eventually it becomes pinned to the lattice 
with a decreased but finite amplitude. However, if (i is 
such that there are sliding velocities in the system, and 
if the soliton starts its motion with the velocity higher 
than some of these then, although it will finally become 
pinned to the lattice, its deceleration will be interrupted 
by long periods of metastable sliding motion at these iso- 
lated velocities. 



B. Concluding remarks 

It is interesting to tie up the above results with our 
previous work on exceptional discretisations of the (j) 4 
model [36]. In that case we discovered that for some 
exceptional models (in which the stationary soliton pos- 
sesses an effective translational symmetry), sliding ve- 
locities, at which the radiation disappears, do exist. The 
fact that all these models involve a complicated nonlo- 
cal discretisation of the nonlinearity leads one to wonder 
whether the nonlincarity has to be discrctised nonlocally 
in order for sliding solitons to exist. This question is 
answered — negatively — by our present work which gives 
a counterexample of a simple, local and physically moti- 
vated nonlinearity supporting sliding solitons. 

A remaining open issue is whether exceptionality of the 
system is a prerequisite for the existence of sliding soli- 
tary waves. Indeed, we have yet to encounter a nonex- 
ceptional discrete system permitting sliding motion. 

We conclude this section by placing our results in the 
context of earlier studies of discrete solitons. 

Moving solitons in the cubic DNLS equation have pre- 
viously been studied by Ablowitz, Musslimani, and Bion- 
dini [13] (among others), using a numerical technique 
based on discrete Fourier transforms as well as pertur- 
bation expansions for small velocities. As we have done, 
these authors suggest that sliding (radiationlcss) solitons 
may not exist. They also observe, with the aid of numer- 
ical simulations, that strongly localised pulses are pinned 
quickly to the lattice, while broader ones are more mobile 
- results which are in agreement with ours. 

Duncan, Eilbcck, Feddersen and Wattis [7, 8] have 
studied the bifurcation of periodic travelling waves from 
constant solutions in the cubic DNLS equation, using nu- 
merical path-following techniques. They find that the 
paths terminate when the soliton's amplitude reaches a 
certain limit value, and in the light of our results it seems 
likely that this is the point at which the radiation be- 
comes large enough to have an effect on the numerics. 
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Solitary waves which decay to constant values at the 
spatial infinities, despite the fact that the generic asymp- 
totic behaviour in the underlying model is oscillatory, are 
commonly referred to as embedded solitons. An exam- 
ple is given by the sliding solitons reported in this paper 
as well as the sliding kinks of [36]; both types of soli- 
tary waves propagate without exciting the resonant back- 
ground oscillations. Embedded solitons have been stud- 
ied for some time in continuous systems [52], but their 
history in lattice equations is younger. Recently, Mal- 
omed, Gonzalez-Perez-Sandi, Fujioka, Espinosa-Ceron 
and Rodriguez have considered certain lattice equations 
with next-to-nearest-neighbour couplings and shown (by 
means of explicit solutions) that both stationary [53] and 
moving [54] embedded solitons exist. Stationary embed- 
ded solitons in discrete waveguide arrays have also been 
analysed by Yagasaki, Champneys and Malomcd [55]. 

Finally, while preparing the revised version of this pa- 
per, we learnt that Melvin, Champneys, Kevrekidis, and 
Cuevas have obtained results very similar to ours. Us- 
ing a combination of intuitive arguments and numerical 
computations, they have found sliding solitons for certain 
values of the parameters of a saturable DNLS model [56]. 
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FIG. 7: Deformation of the integration contour F p to F' p . 



plane; lies in its first quadrant (Rep > 0, Imp > 0); is 
described as a graph of a single- valued function Imp = 
/(Rep) (i.e. never turns back on itself), and is concave-up 
everywhere: 
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We also assume that the function G(p) is analytic in the 
region between the contour 7 and the imaginary axis. 
We begin by writing the left-hand side of (Al) as 



-(p+p"> z F{p)G(p')dp'dp 



- rz G{r-p)dr 



dp, 



(A2) 



APPENDIX A: CONVOLUTION THEOREM FOR 
THE LAPLACE TRANSFORM IN THE 
COMPLEX PLANE 

This appendix deals with the convolution property of 
the modified Laplace transform of the form (37), where 
the integration is along an infinite contour in the complex 
plane rather than the positive real axis. In the context 
of asymptotics beyond all orders, this transform was pi- 
oneered by Grimshaw and Joshi [46]. Since no proof of 
the convolution result (39) is available in literature, and 
since it requires a nontrivial property (concavity) of the 
integration contour, we produce such a proof here. 

We wish to show that 



F{p)dp / e- p ' z G(p')dp' 



F{pi)G(p - pi)dpi 



dp. 



(Al) 



where f£ stands for an integral along the curve 7 from the 
origin to the point p on that curve. We assume that the 
curve 7 extends from the origin to infinity on the complex 



where r = p' +p. The curve T p is the path traced out by 
r as p' traces out the curve 7, for a given p (which also 
lies on 7) - this is depicted in fig. 7. The path T p is the 
same as the path 7 but translated from the origin to the 
point p. 

For each given p we deform the path T p so that it 
now lies on 7, still starting at the point p. We call this 
deformed path r' . The point r — p = p', which lay on 
the path 7 before the deformation, will now move inside 
the region bounded by 7 (i.e. the region to the left of 
7), since all points on T' p lie inside the region bounded 
by r p . This follows from the fact that the path T p is 
concave-up. The point r — p will, however, stay to the 
right of the imaginary axis, since after the deformation all 
points r on T' lie to the right of the point p. This follows 
from the fact that the curve 7 never turns back on itself. 
Since G(p') is analytic for all p between the imaginary 
axis and 7, the value of the integral (A2) will not be 
affected by the deformation. [In the particular case of 
the functions U{p) and W{p) considered in section HIE, 
we have deliberately chosen the contours of integration 
so that the analyticity condition is satisfied.] 

After the deformation, both integrals in (A2) follow 
the path 7, with the inner one starting at the point p. 
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By parameterising the path 7 it is now straightforward 
to change the order of integration, and end up with the 
desired result (Al). 

The argument holds, with the obvious modifications, 
for paths in the second quadrant. 

APPENDIX B: PROOF THAT SINGULARITIES 
DO NOT ACCUMULATE TO THE IMAGINARY 
AXIS 

The roots of eq. (40) with the top and bottom signs give 
singularities of U(p) and W{p), respectively. Our aim 
here is to show that the integration contours j u and 7 S 
in section III E can be chosen to lie above all the complex 
singularities - i.e., that singularities of U (p) and W(p) do 
not accumulate to the imaginary axis. 

The roots of eq. (40) are zeros of the functions F(±p), 
where 

F(p) = coshp — 1 — itanfc(sinhp — p). (Bl) 

The imaginary zeros of F(p) are at p — —iq n , where q n 
are the real roots of eq. (30) . We let N denote the number 
of these imaginary zeros; for k 7^ 0, N is finite. We recall 
that all q n are positive; hence all N imaginary zeros of 
F(p) are on the negative imaginary axis, and all TV zeros 
of F(—p) on the positive imaginary axis. 

Let k and q be the real and imaginary part of p: p = 
K + iq. We let T> denote the rectangular region in the 
complex-p plane bounded by the vertical lines K = ±e 
and horizontal lines q = e 1 ^ 2 and q = Q, where e is small 
and Q is large enough for the region to contain all N 
zeros of F(—p) on the positive imaginary axis. By the 
argument principle, the total number of (complex) zeros 
of the function F(—p) in the region T> is given by (2tt)~ 1 



times the variation of its argument (i.e., its phase) along 
the boundary of T>. In a similar way, one may count the 
number of zeros of the function F(p) in T>. 
We have 

„. . sinliKsino =F tanfc(sinh/t coso — n) 

tan arg F(±p) = y — 7 '— . 

cosh « cos q — 1 ± tan A; (cosh k sin q — q) 

Points on the vertical line k = e satisfy 

tanargF(±p) ps —£-7- In II — cosg ± tanfefg — sing)!. 
dq 

(B2) 

Consider, first, the case of the bottom sign in eq. (B2). 
The expression between the bars in (B2) crosses through 
zero N times as the line k = e is traced out. Each time 
zero is crossed, the logarithmic derivative in (B2) jumps 
from —00 to +00, and the argument of F{—p) increases 
by 7r as we move from one crossing to the next one. The 
net increase of the argument as the line k = e is traversed 
from its bottom to the top, is Nir. In the case of the top 
sign in (B2), on the other hand, the expression in | . . . | 
never crosses through zero and hence the total increment 
of the argument of F(p) is zero. 

As we move along the vertical line n = —s from top to 
bottom, the argument of F(—p) increases by another Nir 
while the argument of F(p) does not acquire any incre- 
ment. Since there are no zero crossings on the horizontal 
segments, the net change of the argument of F(—p) along 
the boundary of V is 27V tt while the total increment of 
arg F(p) is zero. Therefore F(—p) has only N zeros in the 
region T> (and they all lie on the imaginary axis) while 
the function F(p) has no zeros in T>. This implies that 
complex singularities of U(p) and W(p) cannot accumu- 
late to the positive imaginary axis. 
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